pvs[,][nonzeros, i] <- out$pv
pvs[,][-nonzeros, i] <- .5
# if(i %% 100 == 0){cat(i, "; \t")}
}
write.table(inci_cs, file.path(outdir, paste0("X",ix,"_inci.txt")), col.names=F, row.names=F)
write.table(sig_ci, file.path(outdir, paste0("X",ix,"_sigci.txt")), col.names=F, row.names=F)
write.table(pvs, file.path(outdir, paste0("X",ix,"_pvs.txt")), col.names=F, row.names=F)
}
p = 100
outdir <- paste0("checksi_n",n,"_p",p)
dir.create(outdir, showWarnings = F)
ix = 1
set.seed(ix)
X <- matrix(rnorm(n*p),n,p)
beta0 <- c(5:1,rep(0,p-5))
Xb <- X %*% beta0
lambda <- .8/n
inci_cs <- pvs <- sig_ci <- matrix(0, p, ne)   # fraction of times true coefficient in CI
for(i in 1:ne){
eps <- rnorm(n)*1
Y <- Xb + eps
Y <- 1*(Y>mean(Y))
# Fit lasso choosing lambda by CV and min mean AUC, minus 1 to deal with intercept
# cv <- cv.glmnet(x=X[,-1], y=Y, family="binomial",type.measure = "auc")
cv <- glmnet(x=X, y=Y, family="binomial")
beta <- coef(cv, s=lambda, exact=TRUE, x=X, y=Y, family="binomial")   # fixed lambda choice
# lambda <- cv$lambda.min   # lambda chosen
nonzeros <- which(beta[-1] != 0)
out <- fixedLassoInf(x=X, y=Y, beta=beta, lambda=lambda*n, family="binomial", alpha=.05)
inci_cs[,][nonzeros, i] <- 1*(out$ci[,1] < beta0[nonzeros] & out$ci[,2] > beta0[nonzeros])  # true coeff in CIs?
inci_cs[,][-nonzeros, i] <- 1*(beta0[-nonzeros]==0)  # true coeff estimated zero and is actually zero
sig_ci[,][nonzeros, i] <- 1*!(out$ci[,1] < 0 & out$ci[,2] > 0)  # signif nonzero?
sig_ci[,][-nonzeros, i] <- 0  # 1*(beta0[-nonzeros]==0)  # true coeff estimated zero and is actually zero
pvs[,][nonzeros, i] <- out$pv
pvs[,][-nonzeros, i] <- .5
# if(i %% 100 == 0){cat(i, "; \t")}
}
dim(X)
length(beta)
nonzeros <- which(beta[-1] != 0)
out <- fixedLassoInf(x=X, y=Y, beta=beta, lambda=lambda*n, family="binomial", alpha=.05)
p = 10
ix = 1
warnings()
warnings()
assign("last.warning", NULL, envir = baseenv())
warnings()
outdir <- paste0("checksi_n",n,"_p",p)
dir.create(outdir, showWarnings = F)
set.seed(ix)
X <- matrix(rnorm(n*p),n,p)
beta0 <- c(5:1,rep(0,p-5))
Xb <- X %*% beta0
lambda <- .8/n
inci_cs <- pvs <- sig_ci <- matrix(0, p, ne)   # fraction of times true coefficient in CI
for(i in 1:ne){
eps <- rnorm(n)*1
Y <- Xb + eps
Y <- 1*(Y>mean(Y))
# Fit lasso choosing lambda by CV and min mean AUC, minus 1 to deal with intercept
# cv <- cv.glmnet(x=X[,-1], y=Y, family="binomial",type.measure = "auc")
cv <- glmnet(x=X, y=Y, family="binomial")
beta <- coef(cv, s=lambda, exact=TRUE, x=X, y=Y, family="binomial")   # fixed lambda choice
# lambda <- cv$lambda.min   # lambda chosen
nonzeros <- which(beta[-1] != 0)
out <- fixedLassoInf(x=X, y=Y, beta=beta, lambda=lambda*n, family="binomial", alpha=.05)
inci_cs[,][nonzeros, i] <- 1*(out$ci[,1] < beta0[nonzeros] & out$ci[,2] > beta0[nonzeros])  # true coeff in CIs?
inci_cs[,][-nonzeros, i] <- 1*(beta0[-nonzeros]==0)  # true coeff estimated zero and is actually zero
sig_ci[,][nonzeros, i] <- 1*!(out$ci[,1] < 0 & out$ci[,2] > 0)  # signif nonzero?
sig_ci[,][-nonzeros, i] <- 0  # 1*(beta0[-nonzeros]==0)  # true coeff estimated zero and is actually zero
pvs[,][nonzeros, i] <- out$pv
pvs[,][-nonzeros, i] <- .5
# if(i %% 100 == 0){cat(i, "; \t")}
}
warnings()
n
p
n = 200
set.seed(ix)
X <- matrix(rnorm(n*p),n,p)
beta0 <- c(5:1,rep(0,p-5))
Xb <- X %*% beta0
lambda <- .8/n
inci_cs <- pvs <- sig_ci <- matrix(0, p, ne)   # fraction of times true coefficient in CI
assign("last.warning", NULL, envir = baseenv())
for(i in 1:ne){
eps <- rnorm(n)*1
Y <- Xb + eps
Y <- 1*(Y>mean(Y))
# Fit lasso choosing lambda by CV and min mean AUC, minus 1 to deal with intercept
# cv <- cv.glmnet(x=X[,-1], y=Y, family="binomial",type.measure = "auc")
cv <- glmnet(x=X, y=Y, family="binomial")
beta <- coef(cv, s=lambda, exact=TRUE, x=X, y=Y, family="binomial")   # fixed lambda choice
# lambda <- cv$lambda.min   # lambda chosen
nonzeros <- which(beta[-1] != 0)
out <- fixedLassoInf(x=X, y=Y, beta=beta, lambda=lambda*n, family="binomial", alpha=.05)
inci_cs[,][nonzeros, i] <- 1*(out$ci[,1] < beta0[nonzeros] & out$ci[,2] > beta0[nonzeros])  # true coeff in CIs?
inci_cs[,][-nonzeros, i] <- 1*(beta0[-nonzeros]==0)  # true coeff estimated zero and is actually zero
sig_ci[,][nonzeros, i] <- 1*!(out$ci[,1] < 0 & out$ci[,2] > 0)  # signif nonzero?
sig_ci[,][-nonzeros, i] <- 0  # 1*(beta0[-nonzeros]==0)  # true coeff estimated zero and is actually zero
pvs[,][nonzeros, i] <- out$pv
pvs[,][-nonzeros, i] <- .5
# if(i %% 100 == 0){cat(i, "; \t")}
}
warnings()
mean(Y)
hist(Y)
beta <- coef(cv, s=lambda, exact=TRUE, x=X, y=Y, family="binomial")   # fixed lambda choice
beta
n = 100
p
# Test selective inference coverages with systematic simulations
# 4/13/17
#
#### Run this block to test original package
remove.packages("selectiveInference", lib= "~/Documents/BiTEN/selectiveInferenceFM")
install.packages("~/Desktop/selectiveInference", repos = NULL, type="source", lib = "~/Desktop/selectiveInference")
library("selectiveInference", lib.loc = "~/Desktop/selectiveInference")
####
#### Run this block to test Frank's modifications
# remove.packages("selectiveInference", lib= "~/Documents/BiTEN/selectiveInferenceFM")
# install.packages("~/Documents/BiTEN/selectiveInferenceFM", repos = NULL, type="source", lib = "~/Documents/BiTEN/selectiveInferenceFM")
# library("selectiveInference", lib.loc = "~/Documents/BiTEN/selectiveInferenceFM")
####
library("doMC")
library("foreach")
rm(list=ls())
setwd("~/Documents/BiTEN/selectiveInferenceFM")
#### Inputs
nx <- 100
ne <- 1000   # number of errors
sigma <- 1  # error sd
ns <- c(100,1000)
ps <- c(10,50,100,500)
ncores <- 4
####
n = ns[1]
p = ps[3]
outdir <- paste0("checksi_n",n,"_p",p)
dir.create(outdir, showWarnings = F)
ix = 1
set.seed(ix)
X <- matrix(rnorm(n*p),n,p)
beta0 <- c(5:1,rep(0,p-5))
Xb <- X %*% beta0
lambda <- .8/n
inci_cs <- pvs <- sig_ci <- matrix(0, p, ne)   # fraction of times true coefficient in CI
i = 1
eps <- rnorm(n)*1
Y <- Xb + eps
Y <- 1*(Y>mean(Y))
cv <- glmnet(x=X, y=Y, family="binomial")
beta <- coef(cv, s=lambda, exact=TRUE, x=X, y=Y, family="binomial")   # fixed lambda choice
beta
?fixedLassoInf
# Test selective inference coverages
# 4/13/17
#
#### Run this block to test original package
remove.packages("selectiveInference", lib= "~/Documents/BiTEN/selectiveInferenceFM")
install.packages("~/Desktop/selectiveInference", repos = NULL, type="source", lib = "~/Desktop/selectiveInference")
library("selectiveInference", lib.loc = "~/Desktop/selectiveInference")
####
#### Run this block to test Frank's modifications
# remove.packages("selectiveInference", lib= "~/Documents/BiTEN/selectiveInferenceFM")
# install.packages("~/Documents/BiTEN/selectiveInferenceFM", repos = NULL, type="source", lib = "~/Documents/BiTEN/selectiveInferenceFM")
# library("selectiveInference", lib.loc = "~/Documents/BiTEN/selectiveInferenceFM")
####
rm(list=ls())
setwd("~/Documents/BiTEN/selectiveInferenceFM")
#### Inputs
n <- 1000
p <- 50
seed <- 1
ne <- 2e3   # number of errors
mu <- 1     # X mean
sigma <- 2  # X sd
rho <- .5   # correlation in X
####
#### Build X and beta
set.seed(seed)
X <- matrix(NA, n, p-1)
X[,1] <- rnorm(n, mu, sigma)
rho <- .5  # autocorrelated
for(j in 2:ncol(X)){
X[,j] <- rho*X[,j-1] + (1-rho)*rnorm(n, mu, sigma)
}
X <- cbind(1, X)   # with intercept
beta0 <- rep(0,p)
beta0[1] <- -75
beta0[2:6] <- 10
beta0[7:11] <- 1
beta0[12:21] <- .1
Xb <- X %*% beta0
####
set.seed(43)
n <- 1000
p <- 10
X <- matrix(rnorm(n*p),n,p)
# beta0 <- c(5:1,rep(0,p-5))
beta0 <- rep(1,p)
Xb <- X %*% beta0
lambda <- .8/n
inci_cs3 <- pvs3 <- sig_ci3 <- matrix(0, p, ne)   # fraction of times true coefficient in CI
for(i in 1:ne){
eps <- rnorm(n)*1
Y <- Xb + eps
Y <- 1*(Y>mean(Y))
# Fit lasso choosing lambda by CV and min mean AUC, minus 1 to deal with intercept
# cv <- cv.glmnet(x=X[,-1], y=Y, family="binomial",type.measure = "auc")
cv <- glmnet(x=X, y=Y, family="binomial")
beta <- coef(cv, s=lambda, exact=TRUE)   # fixed lambda choice
# lambda <- cv$lambda.min   # lambda chosen
nonzeros <- which(beta[-1] != 0)
out <- fixedLassoInf(x=X, y=Y, beta=beta, lambda=lambda*n, family="binomial", alpha=.05)
inci_cs3[,][nonzeros, i] <- 1*(out$ci[,1] < beta0[nonzeros] & out$ci[,2] > beta0[nonzeros])  # true coeff in CIs?
inci_cs3[,][-nonzeros, i] <- 1*(beta0[-nonzeros]==0)  # true coeff estimated zero and is actually zero
sig_ci3[,][nonzeros, i] <- 1*!(out$ci[,1] < 0 & out$ci[,2] > 0)  # signif nonzero?
sig_ci3[,][-nonzeros, i] <- 0  # 1*(beta0[-nonzeros]==0)  # true coeff estimated zero and is actually zero
pvs3[,][nonzeros, i] <- out$pv
pvs3[,][-nonzeros, i] <- .5
if(i %% 100 == 0){cat(i, "; \t")}
}
####
apply(inci_cs3,1,mean, na.rm=T)
apply(sig_ci3,1,mean, na.rm=T)
apply(pvs3,1,function(z) mean(z < .05, na.rm=T))
beta
source('~/Documents/CSU/17S/STAT675/project/function_file.R')
apply(theta[5000:10000,1:10], 2, mean)
plot(density(theta[5000:10000, 1]))
plot(density(theta[5000:10000, 2]))
plot((theta[5000:10000, 2]), type="l")
acf(theta[5000:10000, 2])
plot(theta[5000:10000, "sigma2"], type="l")
acf(theta[5000:10000, "sigma2"])
plot(theta[5000:10000, "sigma2"], type="l")
acf(theta[5000:10000, "sigma2"])
apply(theta[5000:10000,], 2, mean)
plot(theta[,"sigma2"]*theta[,"tau"]^2, type="l")
mean(theta[5000:10000,"sigma2"]*theta[5000:10000,"tau"]^2, type="l")
source('~/Documents/CSU/17S/STAT675/project/function_file.R')
a1
b1
mean(theta[5000:10000,"sigma2"]*theta[5000:10000,"tau"]^2, type="l")
apply(theta[5000:10000,], 2, mean)
plot(theta[5000:10000, "sigma2"], type="l")
acf(theta[5000:10000, "sigma2"])
acf(theta[5000:10000, 2])
plot((theta[5000:10000, 2]), type="l")
plot(density(theta[5000:10000, 2]))
plot(density(theta[5000:10000, 1]))
1+1/Svec
?rinvgamma
a1 <- 10
b1 <- .1
z <- rinvgamma(1000, shape=a1, scale=b1)
hist(z)
plot(density(z))
mean(z)
a1 <- 10
b1 <- .1
z <- rinvgamma(1000, shape=a1, scale=b1)
hist(z)
plot(density(z))
mean(z)
a1 <- 10
b1 <- 10
z <- rinvgamma(1000, shape=a1, scale=b1)
hist(z)
plot(density(z))
mean(z)
a1 <- 10
b1 <- .1
z <- rinvgamma(1000, shape=a1, scale=b1)
hist(z)
plot(density(z))
mean(z)
a1 <- 10
b1 <- .1
z <- rinvgamma(10000, shape=a1, scale=b1)
hist(z)
plot(density(z))
mean(z)
.1/9
z <- rinvgamma(10000, shape=a1, rate=b1)
mean(z)
1/9/.1
a1 <- 10
b1 <- 10
z <- rinvgamma(10000, shape=a1, rate=b1)
hist(z)
plot(density(z))
mean(z)
10/9
a1 <- 1.1
b1 <- 1.1
z <- rinvgamma(10000, shape=a1, rate=b1)
hist(z)
plot(density(z))
mean(z)
var(z)
a1 <- 1.1
b1 <- 1
z <- rinvgamma(10000, shape=a1, rate=b1)
hist(z)
plot(density(z))
mean(z)
var(z)
a1 <- 11
b1 <- 10
z <- rinvgamma(10000, shape=a1, rate=b1)
hist(z)
plot(density(z))
mean(z)
var(z)
10/9
source('~/Documents/CSU/17S/STAT675/project/function_file.R')
m
source('~/Documents/CSU/17S/STAT675/project/function_file.R')
m
phi
psi
randstart
Svec
b2
b1 + b2/2
a1 + a2/2
a1 + (n+p)/2
source('~/Documents/CSU/17S/STAT675/project/function_file.R')
m
b1/(a1-1)
s2
source('~/Documents/CSU/17S/STAT675/project/function_file.R')
s2
source('~/Documents/CSU/17S/STAT675/project/function_file.R')
m
s2
source('~/Documents/CSU/17S/STAT675/project/function_file.R')
m
s2
?rinvgamma
?rgig
source('~/Documents/CSU/17S/STAT675/project/function_file.R')
source('~/Documents/CSU/17S/STAT675/project/function_file.R')
m
chi
s2
psi
source('~/Documents/CSU/17S/STAT675/project/function_file.R')
m
s2
?rinvgauss
source('~/Documents/CSU/17S/STAT675/project/function_file.R')
m
?rgamma
source('~/Documents/CSU/17S/STAT675/project/function_file.R')
source('~/Documents/CSU/17S/STAT675/project/function_file.R')
m
s2
s <- 1/rgamma(1000, shape=a1 + (n+p)/2,  rate=b1 + b2/2)
mean(s)
r <- rinvgamma(n=1,  shape=a1 + (n+p)/2,  rate=b1 + b2/2 )  # or scale?
mean(r)
source('~/Documents/CSU/17S/STAT675/project/function_file.R')
m
source('~/Documents/CSU/17S/STAT675/project/normal_means_model.R')
source('~/Documents/CSU/17S/STAT675/project/normal_means_model.R')
source('~/Documents/CSU/17S/STAT675/project/normal_means_model.R')
source('~/Documents/CSU/17S/STAT675/project/normal_means_model.R')
plot(theta[,"sigma2"], type="l")
plot(theta[,"a"], type="l")
plot(theta[,"tau"], type="l")
plot(theta[,"sigma2"]*theta[,"tau"]^2, type="l")
plot(table(theta[,"a"]))
outdir
head(theta[,"sigma2"])
apply(theta[5000:10000,1:10], 2, mean)
beta0
plot(density(theta[5000:10000, 1]))
plot(density(theta[5000:10000, 2]))
plot((theta[5000:10000, 2]), type="l")
acf(theta[5000:10000, 2])
plot(theta[5000:10000, "sigma2"], type="l")
acf(theta[5000:10000, "sigma2"])
apply(theta[5000:10000,], 2, mean)
mean(theta[5000:10000,"sigma2"]*theta[5000:10000,"tau"]^2, type="l")
cbind(apply(theta[5000:10000,1:n], 2, mean), apply(theta[5000:10000,1:n], 2, quantile, c(.05,.95)))
apply(theta[5000:10000,1:n], 2, quantile, c(.05,.95))
cbind(apply(theta[5000:10000,1:n], 2, mean), t(apply(theta[5000:10000,1:n], 2, quantile, c(.05,.95))))
beta0
signif <- beta0 > betas[,2] & beta0 < betas[,3]
betas <- cbind(apply(theta[5000:10000,1:n], 2, mean), t(apply(theta[5000:10000,1:n], 2, quantile, c(.05,.95))))
signif <- beta0 > betas[,2] & beta0 < betas[,3]
signif
signif <- 1*(beta0 > betas[,2] & beta0 < betas[,3] )
signif
mean(signif)
source('~/Documents/CSU/17S/STAT675/project/normal_means_model.R')
adistr
tau
la <- suppa*(p*log(tau) - p*log(2) + log(prod(phi))) - p*log(gamma(suppa)) - log(1e8)   # scaled by 1e8
suppa
la
adistr <- ((tau/2)^p * prod(phi) )^suppa / gamma(suppa)^p  # distribution over support of a
adistr
tau
prod(phi)
phi
range(phi)
prod(phi)
m
a <- sample(suppa, 1)
phi <- rnorm(p)^2
phi <- phi/sum(phi)
beta <- rnorm(p)
psi <- rnorm(p)^2
s2 <- rnorm(1)^2
tau <- rnorm(1)^2
a
suppa
suppa <- seq(min(1/n, 1/p), 1/2, length.out=K)    # support of a
a <- sample(suppa, 1)
phi <- rnorm(p)^2
phi <- phi/sum(phi)
beta <- rnorm(p)
psi <- rnorm(p)^2
s2 <- rnorm(1)^2
tau <- rnorm(1)^2
phi
prod(phi)
phi <- rep(1/p,p)
phi
Svec <- psi*phi^2*tau^2
s2 <- 1
sn <- sqrt(1/(1+1/Svec))
for(j in 1:n){
beta[j] <- rnorm(n=1, mean=sn[j]*Y[j], sd=sn[j])
}
for(j in 1:n){
psi[j] <- 1/rinvgauss(n=1,  mean= tau*phi[j]/abs(beta[j]),  dispersion = 1 )
# chi <- (beta[j])^2/s2*tau^2*phi[j]^2
# psi[j] <- rgig(n=1,  chi=chi,  lambda=.5, psi=1 )
}
chi <- 2*sum(abs(beta)/phi)
tau <- rgig(n=1,  chi=chi,  lambda=n*a - n,  psi=1)   # psi is rho in paper
temp <- rep(NA, p)
for(j in 1:p){
temp[j] <- rgig(n=1,  chi=2*abs(beta[j]),  lambda=a - 1,  psi=1)   # psi is rho in paper
}
phi <- temp/sum(temp)
phi
prod(phi)
la <- suppa*(p*log(tau) - p*log(2) + log(prod(phi))) - p*log(gamma(suppa)) - log(1e8)   # scaled by 1e8
adistr <- exp(la)/sum(exp(la))
a <- sample(x=suppa,  size=1,  prob=adistr)
a
adistr <- ((tau/2)^p * prod(phi) )^suppa / gamma(suppa)^p  # distribution over support of a
adistr <- adistr/sum(adistr)   # distribution over support of a, normalized
adistr
a
source('~/Documents/CSU/17S/STAT675/project/normal_means_model.R')
m
prob
adsitr
adist
la
adistr
phi
prod(phi)
randstart
source('~/Documents/CSU/17S/STAT675/project/normal_means_model.R')
m
prod(phi)
phi
hist(phi)
sum(log(phi))
source('~/Documents/CSU/17S/STAT675/project/normal_means_model.R')
a1
betas <- cbind(apply(theta[5000:10000,1:n], 2, mean), t(apply(theta[5000:10000,1:n], 2, quantile, c(.05,.95))))
signif <- 1*(beta0 > betas[,2] & beta0 < betas[,3] )
signif
mean(signif)
source('~/Documents/CSU/17S/STAT675/project/normal_means_model.R')
betas <- cbind(apply(theta[5000:10000,1:n], 2, mean), t(apply(theta[5000:10000,1:n], 2, quantile, c(.05,.95))))
signif <- 1*(beta0 > betas[,2] & beta0 < betas[,3] )
signif
mean(signif)
names(theta)
source('~/Documents/CSU/17S/STAT675/project/normal_means_model.R')
warnings()
source('~/Documents/CSU/17S/STAT675/project/normal_means_model.R')
i
j
r
beta0
Y
hist(Y)
n
source('~/Documents/CSU/17S/STAT675/project/normal_means_model.R')
source('~/Documents/CSU/17S/STAT675/project/normal_means_model.R')
source('~/Documents/CSU/17S/STAT675/project/normal_means_model.R')
m
r
i
j
sigs[1,1]
sigs[1,1,]
mean(sigs[1,1,], na.rm=T)
source('~/Documents/CSU/17S/STAT675/project/normal_means_model.R')
source('~/Documents/CSU/17S/STAT675/project/normal_means_model.R')
source('~/Documents/CSU/17S/STAT675/project/normal_means_model.R')
